Divergent Thermal Conductivity in Three-dimensional Nonlinear Lattices 
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Heat conduction in three-dimensional nonlinear lattices is investigated using a particle dynamics 
simulation. The system is a simple three-dimensional extension of the Fermi-Pasta-Ulam /3 (FPU- 
/?) nonlinear lattices, in which the interparticle potential has a biquadratic term together with a 
harmonic term. The system size is L x L x 2L, and the heat is made to flow in the 2L direction with 
using the Nose-Hoover method. Although a linear temperature profile is realized, the ratio of energy 
flux to temperature gradient shows logarithmic divergence with L. The autocorrelation function of 
energy flux C(i) is observed to show power-law decay as j-°- 9S ±°' 25 ) which is slower than the decay 
in conventional momentum-conserving three-dimensional systems (t~ 3//2 ). Similar behavior is also 
observed in the four-dimensional system. 

PACS numbers: 63.70.+h, 72.25.Dp, 44.10.+i 
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Thermal conduction has been one of the main issues of 
statistical mechanics for more than a century. Thermal 
conduction is usually accurately described by the Fourier 
law 



-kVT, 



(1) 



where J is the heat flux, k is the heat conductivity, and 
T is the local temperature. The transport coefficient is 
described by the Green-Kubo formula[3, H Q 



V 



C{t)At, 



(2) 



where C(t) = (J(t) • J(0)) denotes the equilibrium auto- 
correlation function of J(t). 

In our understanding of the conventional long-time 
tails of the autocorrelation function, C(t) ~ t~ d / 2 in a 
d-dimensional system0,[5, 6j, the integral in eq. (01 is ex- 
pected to diverge in one- and two-dimensional systems, 
and to converge in three- and higher-dimensional sys- 
tems. For a finite system with size L, the size dependence 
of the effective transport coefficient k(L), which is de- 
fined as the ratio of energy flux to temperature gradient, 
is obtained by replacing the upper limit of the integral 
in eq. (0 with the time range L /v s (v s denotes a typical 
phonon velocity). Together with the long-time tail be- 
havior, we obtain k(L) = a + bL 1 ~ d l 2 , where a and b are 
constants depending on the system, and L 1 ~ d / 2 should 
be interpreted to be a logarithmically behaving function 
\nL for d = 2. This argument predicts n(L) ~ L° , \nL 
and a + b/L - 5 for d = 1, 2 and 3, respectively. Therefore, 
k(L) is expected to diverge for d — 1,2 and to converge 
for d> 30,11. 

Such size dependence has been verified in two- and 
three-dimensional fluid systems. In the two-dimensional 
hard-disk system, logarithmic divergence has been con- 
firmed. In the three-dimensional case, 1 / \f~L~ convergence 



has been confirmed in the hard-spheres system 0, 0] and 
in the Lennard- Jones system[lfj ) |ll|. 

In this article, we consider nonlinear lattice systems 
with total momentum conservation as a model sys- 
tem of insulated solids. In one-dimensional nonlinear 
lattices^ S El El El El El, the power-law diver- 
gence of n(L) has been confirmed, where the estimated 
value of the exponent is about 0.37. This value is smaller 
than the one expected from the above argument. In the 
two-dimensional case, logarithmic divergence has been 

observed[HE3- 

The behavior of three-dimensional nonlinear lattice 

models has not been clarified, however. Only a Fermi- 
Pasta-Ulam-/? (FPU-/3)-like three-dimensional model has 
been studied, and 1/y/L convergence was observed|]j. 
This model has a natural length in the potential function, 
and a free boundary condition was used in the direction 
of heat flow. As a result, the system does not have a 
crystalline structure in the steady state. 

The purpose of this study is to determine whether 
three-dimensional nonlinear lattices with momentum 
conserving interaction show 1/y/L convergence. 

Our model is described by the Hamiltonian 
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It is a simple extension of the FPU-/3[19( chain to higher 
dimensional lattices. Here, Pi and j-j denote the mo- 
mentum and displacement of a particle on lattice point i 
respectively, and are three-dimensional vectors. All the 
particles have the unit mass. The summation over 
denotes the nearest-neighbor lattice points. 

A cubic lattice with a size of L x x L y x L z is consid- 
ered. A periodic boundary condition is used in the L x 
and L y directions. A fixed boundary condition is used 
in the L z direction; that is, the particles at both ends 



2 



in the L z direction are coupled to rigid walls through 
the same interaction potential as that between nearest- 
neighbor particles. Furthermore, the particles at both 
ends in the L z direction undergo temperature control by 
the Nose- Hoover method[20|. The temperature at one 
end is denoted by Tl and the other end by Tr; therefore, 
energy flows from the Tl end to the Tr end along the 
L z direction after the system reaches a steady state. In 
summary, the equations of motion are 



Ti = Pi (4a) 

— — — (in the bulk) 

P. { (4b) 
-— dpi (at both ends), 

where Ci denotes the Nose-Hoover thermostat variables, 
which obey 



Q \3k B T 



(5) 



T denotes the target temperatures, Tl and Tji, and 
Q denotes the coupling parameter between Q and pi. 
Boltzmann's constant ks is taken as 1.0. In the following, 
we fix the parameters as 



T L = 20.0, T R = 10.0, and Q = 1.0. 



(0) 



In this temperature region, the nonlinear term in the 
Hamiltonian in eq. has the same order of contribu- 
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FIG. 1: Temperature profile of three-dimensional FPU-/3 
lattice. The size of the system is taken as L x L x 2L. 
Each sequence represents the result for different sizes 2L = 
24, 48, and 96 from top to bottom. The horizontal axis shows 
the z-coordinate rescaled by the system size 2L, and the verti- 
cal axis shows the local temperature averaged over the cross- 
sectional cut in the xy-plane. The 3cr error is within the marks 
for each point. 
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FIG. 2: Size dependence of thermal conductivity in three- 
dimensional FPU-/3 lattices with system size of L x L x 2L 
given in log scale. Logarithmic divergence is observed, as 
shown by the fitted line, (2.047 ± 0.004) x log 2L. 



tion as the linear term. Thus, the temperature is suffi- 
ciently high for the dynamical evolution to reproduce the 
thermal state. 

First, we study the system of size L x L x 2L to avoid 
the effects of anisotropy and dimensional crossover. Sim- 
ulations are carried out using system sizes from L = 4 to 
64. 

Simulations start from a state with all displacements 
r, = and with randomly selected momenta pi, so that 
the local kinetic energy profile satisfies a linear temper- 
ature profile from Tl to Tr. From this initial state, the 
system finally reaches a steady state. This initial relax- 
ation process takes about f = 5x 10 4 for L — 64. After- 
wards, we sample the local temperature T(i z ), where the 



particles are labeled in order as (i x ,i y 



This is given 



by the average of the local temperature T(i) of each par- 
ticle i in the sectional plane of z — i z , with T(i) given by 
the long-time average of the kinetic energy: 



m = <p 4 2 /3). 



(7) 



The temperature profile is shown in Fig. ^ A typical 
simulation time is about 1.0 x 10 5 for L = 64. h There- 
fore, the total number of simulation steps per particle for 
one sample is about 1.0 x 10 7 which takes slightly more 
than 1 month using a single core of a 2.2 GHz Opteron 
processor for the L = 64 system. Five to eight samples 
are accumulated for the results. The temperature profile 
T(z) becomes linear, which shows that thermalization is 
sufficient in this state. 

We estimate the local heat flux at site i with the energy 
flow using 



1 



(ri+i,+r i )--^ r V(r i+u -r i ), (8) 
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FIG. 3: Autocorrelation function C(t) = (J z (t)J z (0)} in 
three-dimensional FPU-/?-lattice with system size of 32 x 
32 x 64. The temperature is set to T — 15. The peri- 
odic boundary condition is applied, and the power-law ex- 
ponent of the long-time tail behavior is shown by the fitted 



line (C(t) ~ t 



-0.98±0.25 



where i + l z denotes the nearest-neighbor site of the zth 
particle in the z-direction. Heat flux per particle, J z , is 
estimated using 



Jz 




(9) 



where (•) denotes the time average after the simulation 
reaches to a steady state and N denotes the total num- 
ber of particles. Then, the thermal conductivity k(L) is 
estimated using 



"'" = ■'■/(§) 



(10) 



Figure |21 shows the estimated values of this n(L). A log- 
arithmic divergence is clearly observed for the systems 
up to 64 x 64 x 128. Divergence of k(L) with L means 
that the system does not have finite conductivity at the 
macroscopic limit. 

This logarithmic divergence of k(L) is consistent with 
the long-time behavior of the autocorrelation function, 
C(t) = (J z (t)J z (0)) in the equilibrium state. Figure 
shows the estimated values of this C (t) obtained by 
microcanonical simulation of the same nonlinear lat- 
tice without temperature control. The system size is 
32 x 32 x 64, and the total energy in the system is adjusted 
to the internal energy expectation value at temperature 
T = 15. A periodic boundary condition is used in all di- 
rections in this simulation. In Fig. 3, it is observed that 
C(t) decays asymptotically as i/£°- 98 ± - 25 m the long- 
time limit, which suggests the logarithmic divergence of 
k(L), as we observed in our nonequilibrium simulation. 



FIG. 4: Size dependence of thermal conductivity in quasi-one- 
dimensional FPU-/3 lattices. Circle, triangle, and square plots 
respectively represent the thermal conductivities for system 
sizes of 3 x 3 x L, 4 x 4 x L, and 8 x 8 x L. 



This is evidence that the divergence of k(L) is simply due 
to the bulk property of the system, and it is not caused 
by a boundary effect or by the temperature control. 

Secondly, the thermal conductivity of quasi-one- 
dimensional systems is studied. For fixed L x and L y , 
the L z dependence of the thermal conductivity, k(L z ), is 
estimated. Details of the computer simulation and the 
model parameters arc the same as above. The systems 
with L x = L y =3, 4, and 8 are simulated for L z = 8 
to 1024. The results are shown in Fig. 0] Estimated 
values of thermal conductivity are consistent with the 
bulk values of the system with the same L z as shown in 
Fig. 2, for up to L z ~ 128, and the longer systems have 
larger k and show power-law divergence, which was con- 
firmed in various one-dimensional systems. The crossover 
length from the logarithmic divergence to the power-law 
increases with thickness (L x = Ly). The length is about 
256 or 512, even for the thinnest L x = L v = 3 system. 
Therefore, the crossover length is about 100 times the 
cross-section size. This result may explain the reason why 
the convergence of thermal conductivity was confirmed 
in previous studies that used quasi-one-dimensional sys- 
tems. 

So far, we have observed that conductivity is logarith- 
mically divergent not only in two-dimensional, but also in 
three-dimensional FPU-/3 lattices. We show here that the 
conductivity of the four-dimensional FPU-/? lattice also 
has the tendency of logarithmic divergence. We prepare 
the hypercubic lattice with a system size of L x L x L x 2L 
in four-dimensional space, and simulate the system using 
eqs. ©, l|4a|l. and (|4bjl . and boundary conditions similar 
to those used in the three-dimensional model. The same 
parameters are used as in the three-dimensional case ex- 
cept for the temperatures, which are set to be Tl = 15.0 
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and Tr = 7.5. The calculated thermal conductivity is 
shown in Fig. [5] It is notable that the tendency of di- 
vergence still appears, as in the three-dimensional case, 
although we could only calculate up to the system size of 
2L ~ 32 because larger systems require too much com- 
putational load. The result seems to indicate that no 
dimensionality effect exists in the behavior of the con- 
ductivity at the isotropic thermodynamic limit when the 
dimension is higher than 2. 

The behavior of the model with a natural length and 
fixed boundary conditions has not yet been investigated, 
although the convergence of k{L) has been observed in 
the free-boundary-applied FPU-/3 lattice with a natural 
length.0. We simulated such systems with a size of 5 x 
5 x L. We set the natural length as Iq — 100.0, and the 
boundary condition in the x and y directions as periodic, 
and in z direction, the particles at both ends are linked 
to missing atom which cannot move. The missing atoms 
on the left and the right sides are separated by a distance 
of (L + 1) X Iq. In such a situation, the conductivity n(L) 
is confirmed to show diverging behavior for a system size 
of up to L = 512. 

In this study we investigated heat transport in nonlin- 
ear lattices with momentum conservation as a model of 
insulated solids. The model we used contained no impu- 
rities or randomness in mass or interaction, and only the 
effects of nonlinear interaction and dimensionality con- 
tribute to thermalization. We find that such nonlinear 
lattice systems do not have finite thermal conductivity 
at the thermodynamic limit, even if they are three- and 
four-dimensional systems. Contrary to the standard un- 
derstanding of t~ d / 2 long-time tail behavior, £ _1 behavior 
is widely observed not only in two-dimensional lattices, 
but also in three- and four-dimensional lattices, and may 
also appear in higher-dimensional lattices. Therefore, the 




standard long-time tail argument cannot be applied in 
the case of nonlinear lattices. It is natural, because the 
conventional long-time tail is considered a consequence 
of viscous modes, which do not exist in our model solids. 
Our next challenge is to clarify the kinematical origin 
of the t~ x behavior. Of course, studying systems with 
larger sizes is a topic of future study. 

We have not investigated in detail how the conduc- 
tivity behaves when we increase the values of nonlinear 
strength g and temperature T. However, we confirm a 
similar divergence of k(L) when g is up to 1.0 and T is 
up to approximately 100.0. Detailed investigations on 
the effect of stronger nonlinearity are problems to be in- 
vestigated. 

It is clear that there is still a long way to go before 
we can model normal heat conduction using nonlinear 
lattice models. Or, our results may be an indication of 
a new mechanism of heat flow. The problems concern- 
ing the minimum conditions for normal heat conduction 
and the behavior of heat flow in insulated solids remain 
unsolved. Heat conduction in three-dimensional mod- 
els requires further investigation to realize a conclusive 
model for heat conduction in crystalline insulated solids. 
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FIG. 5: Size dependence of thermal conductivity in four- 
dimensional FPU f3 lattices using semilog scale. The size of 
the system is L x L x L x 2L. 



[1 

[2 

[3] 
[4] 

[5 

[f 

[8] 
[9] 

[io; 
in 

[12 
[13 



* Electronic address: shibaQacolyte.t.u-tokyo. ac.jp 

' Present address: Department of Earth and Space Sci- 
ence, Graduate School of Science, Osaka University, Os- 
aka 560-0043 

* Electronic address: ito@ap.t.u-tokyo. ac.jp 
R. Kubo: J. Phys. Soc. Jan. 12 (1957) 570. 
R. Kubo, M. Yokota and S. Nakajima: J. Phys. Soc. Jpn. 
12 (1957) 1203. 

M. S. Green: J. Chem. Phys. 22 (1954) 398. 
B. J. Alder and T. E. Wainwright: Phys. Rev. A 1 (1970) 
18. 

Y. Pomeau and P. Resibois: Phys. Rep. 19 (1975) 63. 
M. H. Ernst: Physica D 47 (1991) 198. 
T. Shimada, T. Murakami, S. Yukawa, K. Saito and N. 
Ito: J. Phys. Soc. Jpn. 69 (2000) 3150. 
S. Lepri, R. Livi and A. Politi: Phys. Rep. 377 (2003) 1. 
T. Murakami, T. Shimada, S. Yukawa and N. Ito: J. 
Phys. Soc. Jpn. 72 (2003) 1049. 

F. Ogushi, S. Yukawa and N. Ito: J. Phys. Soc. Jpn. 74 

(2005) 827. 

F. Ogushi, S. Yukawa and N. Ito: J. Phys. Soc. Jpn. 75 

(2006) 073001. 

H. Kaburaki and M. Machida: Phys. Lett. A 181 (1993) 
85. 

S. Lepri, R. Livi and A. Politi: Phys. Rev. Lett. 78 (1997) 



■5 



1896. 

[14] T. Hatano: Phys. Rev. E 59 (1999) Rl. 
[15] S. Lepri, R. Livi and A. Politi: Europhys. Lett. 43 (1998) 
271. 

[16] O. Narayan and S. Ramaswamy: Phys. Rev. Lett. 89 
(2002) 200601. 

[17] A. Lippi and R. Livi: J. Stat. Phys. 100 (2000) 1147. 



[18] S. Lepri, R. Livi and A. Politi: Chaos 15 (2005) 015118. 

[19] E. Fermi, J. Pasta and S. Ulam: Los Alamos Report 
No. LA-1940 (1955), later published in E. Fermi, Col- 
lected Papers, E. Segre (ed.), University of Chicago Press 
(1965). 

[20] S. Nose: Prog. Theor. Phys. Suppl. 103 (1991) 1. 



